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ABSTRACT 

X-ray afterglow light curves have been collected for over 400 Swift gamma-ray 
bursts with nearly half of them having X-ray flares superimposed on the regular 
afterglow decay. Evidence suggests that gamma-ray prompt emission and X-ray 
flares share a common origin and that at least some flares can only be explained 
by long-lasting central engine activity. We have developed a shell model code to 
address the question of how X-ray flares are produced within the framework of the 
internal shock model. The shell model creates randomized GRB explosions from 
a central engine with multiple shells and follows those shells as they collide, merge 
and spread, producing prompt emission and X-ray flares. We pay special atten- 
tion to the time history of central engine activity, internal shocks, and observed 
flares, but do not calculate the shock dynamics and radiation processes in detail. 
Using the empirical Ep — Eiso (Amati) relation with an assumed Band function 
spectrum for each collision and an empirical flare temporal profile, we calculate 
the gamma-ray (Swift/BAT band) and X-ray (Swift/XRT band) lightcurves for 
arbitrary central engine activity and compare the model results with the obser- 
vational data. We show that the observed X-ray flare phenomenology can be 
explained within the internal shock model. The number, width and occurring 
time of flares are then used to diagnose the central engine activity, putting con- 
straints on the energy, ejection time, width and number of ejected shells. We 
find that the observed X-ray flare time history generally reflects the time history 
of the central engine, which reactivates multiple times after the prompt emission 
phase with progressively reduced energy. The same shell model predicts an ex- 
ternal shock X-ray afterglow component, which has a shallow decay phase due 
to the initial pile-up of shells onto the blast wave. However, the predicted X-ray 
afterglow is too bright as compared with the observed flux level, unless e^, is as 
low as 10~^. 

Subject headings: gamma-rays: bursts — shock waves 
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1. Introduction 

The study of gamma-ray bursts (GRBs) has been greatly advanced following the launch 
of the Swift Gamma-Ray Explorer (Gehrels et al. 2004) on November 20, 2004. Thanks 
to its rapid slewing capability, multi-wavelength observations of GRB afterglows, usually 
as soon as < 100 seconds after the burst trigger, have been performed regularly for many 
bursts. The X-ray telescope (XRT, Burrows et al. 2005a) aboard this satellite has given 
us unprecedented access to early afterglows in the X-ray band. For most GRBs, XRT has 
observed smoothly decaying power-law or broken power-law afterglows (Nousek et al. 2006; 
O'Brien et al. 2006; Zhang et al. 2007; Liang et al. 2007, 2008; Evans et al. 2009). In about 
half of GRBs, superimposed on the background X-ray afterglow, one or more X-ray flares 
are seen (Burrows et al. 2005b; Romano et al. 2006; Falcone et al. 2006, 2007; Chincarini 
et al. 2007). 

Flares among the bursts share some common properties, which include the following: 

• The morphology of flares is similar: they all show a smooth, rapid rise and a rapid fall 
(Romano et al. 2006). 

• They are superimposed on a background power-law decay afterglow component with 
the slope before the flaring equal to that after the flare (Burrows et al. 2005c). 

• The width of flares is typically narrow, with 5t/t ~ 1/10 on average, where t is the 
emission time of the flare. However, flares can become wider {5t/t becomes larger) 
as the time of their emission increases, and no sharp flares are seen at late times 
(Chincarini et al. 2007; Kocevski et al. 2007). 

• Similar flaring activity has been detected in both types of GRBs: those believed to be 
of massive star origin (Type II, typically long duration) and those believed to be of 
compact star merger origin (Type I, typically short, e.g. in GRB 050724, Barthelmy 
et al. 2005)0. This suggests that the phenomenon is insensitive to the progenitor type 
(Perna et al. 2006). 

• Flares are spectrally harder than the underlying afterglow (Burrows et al. 2005b; 
Romano et al. 2006; Falcone et al. 2006). 

• For a small sample of bright flares whose time- dependent spectral analysis can be 
performed, flares are found to soften as they decay (Burrows et al. 2005b, c; Falcone 
et al. 2006, 2007), reminiscent of GRB prompt emission. 



^For a full discussion of the two physically distinct types of GRBs, see Zhang et al. (2009). 
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• Average X-ray flare luminosity decays with time as a power- law with slope ~ -1.5 for 
a sample of flares from many GRB's. This temporal relationship also seems to hold 
among flares from a single multi- flared burst (Lazzati, Perna & Begelman 2008). 

Flares also vary from burst to burst, from flare to flare, in the following ways: 

• The shape of flares can be best flt with a number of different proflles: gaussian, log- 
normal distributions, exponential or power-law rise or decay with differing slopes for 
both the rising and falling portion of the flare (Chincarini et al. 2007). 

• The number of flares seen per burst varies from modal value of 1, mean value of ~ 2.5 
and maximum value of 8 (Chincarini et al. 2007). 

• The fluence seen in flares can vary from a few percent of the burst fluence to larger 
than the burst itself (e.g. for GRB050502B, Burrows et al. 2005b, c; Falcone et al. 
2006). 

• The emission time of flares varies from perhaps before the slew of XRT (e.g. GRB060607) 
to late flares at 10^ seconds (e.g. GRB050904, Cusumano et al. 2006) to 10^ seconds 
(for the Type I GRB050724, Barthelmy et al. 2005). Some flares are superimposed 
on other flares (e.g. GRB050916, Burrows et al. 2005b,c; Chincarini et al. 2007). 
Although most flares happen early, from 100 to 1000 seconds, the distribution of t tails 
off to 10^ seconds (Chincarini et al. 2007). 

A few properties of X-ray ffares suggest that they are connected to internal emission 
processes due to a restarting of the central engine at late times (Burrows et al. 2005b; 
Falcone et al. 2006; Romano et al. 2006; Zhang 2007). The arguments in favor of such 
a "late internal" model (in contrast to the external shock model) are the following: First, 
the rapid rise and fall of flares with 5t/t ~ 1/10 strongly disfavors external shock models 
that involve a large angular area over which radiation would be emitted (Zhang et al. 2006, 
see also loka et al. 2005; Fan & Wei 2005; Lazzati & Perna 2007). Second, the connected, 
underlying continuum which has the same slope both before and after the flare suggests 
that flares are not related to the canonical afterglow and are instead superimposed on top of 
this regular afterglow decay (Chincarini et al. 2007). Third, given the same observed X-ray 
flare amplitude, the internal model requires a much smaller energy budget than the external 
shock model (Zhang et al. 2006). While the internal model can produce signiflcant X-ray 
flares with energy much less than that during the prompt emission, the external shock model 
requires an energy budget at least comparable to that of the initial blast wave in order to 
make a noticeable change in flux (Zhang & Meszaros, 2002). Internal models are much more 
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"economical" as far as energy budget is concerned. Next, Liang et al. (2006) have shown that 
as long as the central engine clock is reset to zero at the beginning of each flaring episode, 
the curvature effect can naturally explain the spectral index and the temporal decay index 
of the decay phase of flares. Finally, although there is no correlation between the number of 
prompt emission burst pulses and X-ray flares in any given burst (Chincarini et al. 2007), 
X-ray flares do exhibit the same spectral softening as GRB prompt emission (Burrows et al. 
2005b,c; Falcone et al. 2006, 2007). This suggests that the properties of X-ray flares make 
them likely to be caused by the same mechanism as prompt emission seen in gamma-rays. 

Any X-ray flare model must be able to explain the bulk similarities and differences 
among flares enumerated above. The leading "internal" model is the internal shock model 
(Rees & Meszaros 1994), which invokes internal collisions of shells within an unsteady central 
engine wind. Previous internal shock models have focused on interpreting prompt gamma- 
ray emission properties (e.g. Kobayashi et al. 1997; Panaitescu et al. 1999; Spada et al. 
2000; Guetta et al. 2001). Kobayashi et al. (1997) have shown that internal shocks can 
explain the highly variable profiles seen in GRBs. In their model, shells of matter with 
varying Lorentz factors ejected from the central engine produce about as many collisions 
as the number of shells, and the time of ejection of the shell is highly correlated to the 
time of collision of the shells, indicating that there is essentially no delay between energy 
ejection by central engine emission and the time when that emission is seen. Recently, using 
a model where coUisions are not perfectly inelastic, Li & Waxman (2008) studied the effect 
of residual collisions in optical emission which would be slightly delayed from gamma-ray 
prompt emission. 

Within the context of X-ray flares, the internal shock model has been discussed by a 
number of authors. Zhang et al. (2006) and Fan & Wei (2005) discussed how a late internal 
shock may produce a softer flare than prompt gamma-ray emission. Wu et al. (2005) 
discussed both late internal and external shock models to study X-ray flares and concluded 
that at least some flares have to be produced by late internal shocks. Lazzati & Perna (2007, 
see also Zhang 2007) proved the suggestions of Burrows et al. (2005b) and Zhang et al. 
(2006) that late flares must require late injection of shells and cannot be produced by late 
collisions of shells ejected during the prompt phase. Yu & Dai (2009) studied the shock 
physics and radiation processes of late internal shocks that may be responsible for X-ray 
flares. 

In this paper, we focus on another aspect of the internal shock model for X-ray flares. 
Extending the work of Kobayashi et al. (1997) to concentrate on X-ray flares, we have created 
a GRB fireball model with a central engine that can eject multiple episodes of matter shells 
with any Lorentz factor, thickness and mass distributions. This code is used to address the 
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question of what kind of central engine activities are demanded in order to reproduce the 
observed properties of X-ray flares. The conclusion drawn from this study may be taken 
as the requirements of some central engine models on X-ray flares (e.g. King et al. 2005; 
Perna et al. 2006; Proga & Zhang 2006; Dai et al. 2006; Lee et al. 2009)1. In §2, we first 
examine the physics of two-shell collisions and describe our phenomenological treatment of 
the spectrum and lightcurve of individual colhsions. In §3, we discuss blast wave evolution, 
which sets an outer boundary for the internal collisions that are relevant to X-ray flares. 
We then dedicate our discussion (§4) to multiple collisions for multiple ejection episodes 
and use the observations to diagnose the required central engine activities. Our results are 
summarized in 55 with some discussion. 



Our model represents the GRB central engine as a variable source that ejects many 
randomized shells of matter over some period of time. Shells ejected later with higher 
Lorentz factors catch up with slower shells, colliding and creating "internal shocks" (Rees & 
Meszaros, 1994). As the modeled shells move outward from the central engine, the energy 
released from each collision is calculated, allowing the two colliding shells to combine and 
continue to move outward to collide with other shells. In the model, each shell is given four 
initial parameters: ejection time tej, relativistic Lorentz factor 7, mass m, and thickness A. 

We assume that shell collisions are inelastic, so that a fast shell (f) and a slow shell (s) 
merge to form a merged shell (m). Using conservation of momentum and energy, the Lorentz 
factor of the combined shell can be written 



2. 



Two-Shell Interaction 




(1) 



Each collision then releases an internal energy given by (Kobayashi et al. 1997) 



Eint = (7/ - 1m)mfC^ + (75 - 7m)msC^. 



(2) 



The efficiency of each collision can be defined as 



7] 



(3) 



^ Other internal dissipation models for X-ray flares (e.g. Panaitescu 2008) may not subject to these 
requirements. 
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It has been known that this efficiency is usually low (Kumar 1999; Panaitescu et al. 1999) 
although depending on input parameters, it can vary in a wide range. To test this, we 
simulated 100 shells that are randomly injected during 0-100 s. We allow the mass of each 
shell to be randomly drawn in the range of 10^^ — 10^^ g in log space, and investigate how 
Lorentz factor contrast affects the distribution of the efficiency. We randomly generate shell 
Lorentz factors in log space within a range of (7™™, Imax)- Table 1 shows the mean efficiency 
of 77 and its standard deviation a for varying 'jmax/lmin- It is evident that the mean value 
of 77 rises steadily with ^max/lmin- However, even for a very high Lorentz factor contrast 
Imaxhmin — 1000, the mean efficiency still only reaches a level of ~ 26%. Values of a are 
on the order of 77, which means that very low and very high -q can be expected for extreme 
parameters of the two shells. The right pane of Table 1 shows the efficiency for the same 
simulation, but with shell masses standardized to all be 10^° g. These results are on order 
with the efficiency study with varying shell masses. 

For our model, we do not follow the internal shock physics in detail. This has been 
done by Yu & Dai (2009), who showed that for reasonable parameters, the emission from 
internal shocks can well reproduce the X-ray flare phenomenology. Our goal is to model 
many collisions and to investigate the time history of central engine activity, collision, and 
GRB emission. Like many other previous work in this direction (e.g. Kobayashi et al. 1997; 
Panaitescu et al. 1999), we adopt an empirical approach to calculate the spectrum and 
lightcurve of each coUision event. 

Spectral model. The spectrum of GRBs is typically a smoothly-joint-broken-power-law 
spectrum, or "Band" -function (Band et al. 1993) 

/ ^(lol^rexp(-t) {a-m>E 

Ne{E) = { , xa-B (4) 

1 ^ (^) exp(/? - a) {-^Y {a - (3)Eo < E 

where a and j3 are the power law photon spectral indices below and above the break energy 
Eq in the asymptotic regime, A is a normalization factor. The break energy can be written 



Table 1: Efficiency and standard deviation for energy conversion in the internal shock model. 



'ymax / '~)min ^ 


'^max / 'ymin ^ ^ 


variable masses 


even masses 


1000 26% 26% 
100 14% 17% 
10 5% 7% 
5 0.5% 0.7% 


1000 28% 27% 
100 19% 20% 
10 7% 8% 
5 0.6% 0.9% 
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as Eo — Ep/{2 + a), where Ep is the spectral peak energy in the GRB energy spectrum. 
Since evidence strongly suggests that GRB prompt emission and X-ray flares originate from 
similar physical events (e.g. Burrows et al. 2005c; Chincarini et al. 2007; Falcone et al. 2007; 
Krimm et al. 2007), we assume that the spectrum of X-ray flares is also a Band-function. 

The Band-function parameters of X-ray flares are determined in the following way in our 
calculations: The typical values for the spectral indices are taken as a = —1, and /3 = —2.3. 
For Ep, we calculate it through an empirical relation between the isotropic emission energy 
Eiso and Ep, which is generally valid among GRBs (Amati et al. 2002; Krimm et al. 2009) 
and within a burst (Liang et al. 2004; Ghirlanda et al. 2009). We assume the validity of 
this correlation 



and apply Eiso — E^nt to estimate Ep. This is because electrons in the shock are in the 
"fast-cooling" regime, and lose their energy rapidly. As long as the electron cquipartition 
parameter eg is close to unity, essentially all the internal energy can be radiated away. For 
eg ^ 1, this estimate gives an upper limit on the brightness of X-ray flares. We require 



for the normalization factor A. Then reinserting this constant, we can calculate the energy 
of each pulse ENE{E)dE within each band {Ei, E2), e.g. 15 to 150 keV for BAT, 0.1 to 
10 keV for XRT. 

Temporal model. The total energy released in a particular band for a particular pulse 
is distributed in time throughout the X-ray flare temporal proflle. Different functions of 
GRB pulses and X-ray flares have been adopted in the hterature (Kobayashi et al. 1997; 
Chincarini et al. 2007), but for the purpose of our study, the shape of the flare proflle is not 
crucial. For simplicity, we adopt the following temporal proflle 



where the temporal width St — A/c scales as the physical width of the shell, A. For a shell 
with initial width Aq moving outward from the central engine, spreading will occur after the 
sound wave travels across the shell at Rg ~ 7^Ao). So in general, the shell width can be 
expressed as (Meszaros et al. 1993; Kobayashi et al. 1999) 




(5) 



that the spectrum-integrated bolomctric energy ENE{E)dE equals Eint and then solve 



L{t) = L{tp)e'^(^ 



(6) 




R< R 
R> R 



(7) 



We require that the temporal integral of the proflle L{t)dt equals the internal energy 
in the specifled energy band (BAT or XRT). When shells collide, the duration of the X-ray 
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flare pulse is defined by the width of the faster shell. This is because the reverse shock 
is typically the one that dominates the X-ray flare emission (Yu & Dai 2009). After the 
collision, which occurs at Rcoi, the width of the the combined shell is taken to be Rcoi/lm^ 
i.e. the width in the spreading regime. The width of the shell keeps spreading as -R/t^ 
thereafter. 

In the central engine frame, the shells collide at t^oi — Rcoi/c- The relevant observation 
time i^ 

. _ J. I ^^(^ol — tej) /ON 
r®,coi — tej -\ ; \0) 

where tej and 7 can be taken as the ejection time and Lorentz factor of either of the two 
colliding shells. The peak time of a flare is defined as the observed collision time plus the 
observed shock crossing time, which is roughly estimated as 

tp ^ t(^,col + — ■ (9) 



3. Blast Wave Evolution 

In order to define "internal" collisions, we need to track the location of the "external" 
shock. As the first shell moves outward into the ambient medium, it slows down due to an 
external shock mechanism as it sweeps up the ambient medium. As time goes by, more and 
more trailing shells collide onto the leading decelerating shell. The motion of this decelerating 
ejecta along with the medium collected along the way, known as the "blast wave" , is governed 
by the following differential equations (Chiang & Dermer 1999; Huang et al. 2000): 



= f3c= c, 10 

at 7 



±L = =t^, (u) 

dm M ' ^ ' 



•^Usually this is called observer frame time. However, the central engine and the observer is in the same 
inertial frame (with cosmological time dilation correction). The difference between the two times is due to 
a propagation effect, not Lorentz transformation (Zhang & Meszaros 2004). 
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Here t is the time in the rest frame of the central engine, R is the distance from the central 
engine, p is the density of the ambient medium, 7 is the Lorentz factor of the shell, m is the 
swept-up mass, and M = Mq + 7m is the total mass including internal energy of the blast 
wave, where Mq is the initial mass of the ejecta. As a result, one has another differential 
equation 

dM 

Equation ffTOl) simply states how radius changes as a function of time for an object moving 
with constant velocity. Equation ffTTl) is a statement of conservation of energy and momentum 
across the blast wave (Blandford and McKee 1976). Although these expressions are valid 
only when the blast wave is in the relativistic regime, they are adequate for our calculation, 
since X-ray flares usually happen early before the blast wave enters the trans-relativistic 
regime. The amount of swept up mass in shell of surface area AttR^ is described by equation 
(fT2ll R. The solution of this system is found to be 



7 



/ IGn^pR^ + 247rp(Mo7o)i?3 + 9(Mo7o)2 
y levrV^e + 247rp(Mo7o)i?3 + 9M2 ' ^ ^ 



The leading shell initially expands freely until the momentum of the swept up matter 
is about equal to the initial mass of the shell, Mq. The radius of deceleration, Ra, is found 
from 'Jq^hR^p ^ Mq, where the subscript "0" represents initial values (Rees and Meszaros, 
1992). For typical values of 70 ~ 100, Mq ~ lO^^g and p ~ 2 x 10~^^g cm~^, one has 

■d 



~ 2 X lO^^cm. 



The relativistic blast wave decays as 7 cx R~^^'^, until reaching the Sedov radius, at 
which the rest mass of the ambient medium becomes as large as the rest mass of the blast 
wave, i.e. Rsedov = (3Mo/47rp)^/^. After this point, the expansion enters the non-relativistic 
regime with dimensionless velocity /? oc R~^^^ (Sedov 1969). 

In our model, trailing shells collide amongst themselves and later land onto the blast 
wave as the it slows down, altering the blast wave dynamics. The dynamics of such a collision 
is very complicated, invoking three shocks and several distinct dynamical stages (Zhang & 
Meszaros 2002). For the purpose of this study (tracking the location of the blast wave), 
we adopt the following simple treatment: If a trailing fast shell with mass rrif and Lorentz 



^Tlic treatment here is based on the assumption of an isotropic ejecta. The treatment is valid before 
the "jet break" time, which is usuaUy the case for X-ray flare observations. After the jet break time, the 
dynamics may be altered by sideways expansion of the ejecta. However, numerical simulations suggest that 
such an expansion is not important (Kumar & Granot 2003; Cannizzo et al. 2004; Zhang & MacFadyen 
2009). For simplicity, we take the isotropic assumption throughout the blast wave evolution. 
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factor 7/ collide on to the blast wave with Lorentz factor of 7 before the blast wave has begun 
decelerating, the collision is treated as an internal shock and the merged Lorentz factor are 
calculated as prescribed by Eq.([T]). If the collision of a fast shell onto the blast wave occurs 
after the blast wave has begun to decelerate, we first assume that there is no blast wave, and 
calculate the post-collision product the same way as the internal shock calculations, record 
the new effective initial mass and Lorentz factor, and re-solve the blast wave for the new 
parameters. We then jump the blast wave solution from the old (low) 7 value to the new 
(high) 7-value at the same R. The code then tracks this new solution until next collision 
happens. With each new collision onto the blast wave, we then calculate the new solution 
by changing the effective mass and initial 7 factor. By doing so, the blast wave evolution 
shows several "glitches" with decreasing amplitude, since the ratio between the trailing shell 
energy and the blast wave energy drops with time as the energy of the blast wave grows (see 
Figure [T]) . 

The location R and Lorentz factor 7 of the blast wave are traced at any instant t. This 
information is used to screen out the collisions that are not "internal". Without such a 
blast wave screening, shells injected from the central engine can collide at any location at 
any time. If the relative Lorentz factor of two shells is small, they can in principle collide 
at a much larger radius than the blast wave radius. Should these collisions have happened, 
both shells have already entered the "spreading" regime, so that the width of the shells can 
become very large. This would introduce some "fat" X-ray flares which are not observed. 
The blast wave essentially restricts such collisions, ensuring that they would never happen. 
Instead, both shells collide onto the blast wave and boost up the blast wave energy. The 
inclusion of the blast wave dynamics is therefore essential, screening out many "fat" flares. 
This also explains the lack of very "fat" flares in the X-ray afterglow data. 

In order to reproduce the data, we adopt a typical afterglow template according to the 
data (Nousek et al. 2006; O'Brien et al. 2006; Liang et al. 2007) to denote the underlying 
X-ray afterglow. This is a broken power law with decay indices -1/2 and -1.25, with a break 
time tfe = 6 X lO^s and a break flux Flux(tb) = 1.3 x 10~^^ erg s~^ cm~^. The origin of the 
X-ray afterglow, especially the shallow decay phase, still remains a mystery. A number of 
very different physical mechanisms have been proposed to explain its origin (see Zhang 2007 
for a review). Although the refreshed shock model is mostly discussed (Zhang et al. 2006; 
Nousek et al. 2006; Granot & Kumar 2006), some optical afterglows do not show a similar 
temporal break around the X-ray break time (Panaitescu et al. 2006; Liang et al. 2007) 
suggesting that this model cannot interpret all the X-ray afterglow data. Other ideas/models 
include a central engine powered afterglow (Ghisellini et al. 2007; Kumar et al. 2008), a 
long-lasting reverse-shock-dominated afterglow (Genet et al. 2007; Uhm & Beloborodov 
2007), two-component external shock (Racusin et al. 2008; De Pasquale et al. 2009), dust 
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10l5 iQl6 iQl7 iQl8 iQl9 iq20 

R [cm] 

Fig. 1. — The blast wave shows several glitches as more and more shells pile onto it. Succes- 
sive solutions are shown as dashed lines and the blast wave is shown as a solid line. Jumps 
in the solution become less important as the added energy from a single collision becomes 
less of a fraction of the total energy. 
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scattering (Shao & Dai 2007), up-scattering of blast wave photons by a trailing lepton-rich 
ejecta (Panaitescu 2008) and an emission component prior to the GRB trigger (Yamazaki 
2009; Liang et al. 2009). Since none of the scenarios have been robustly proven, we do not 
demand our model to interpret the power law X-ray afterglow segments self-consistently. 

On the other hand, the expected external shock X-ray afterglow flux can be calculated 
from our model. Applying the standard afterglow model (Sari et al. 1998), we can calculate 
the X-ray lightcurve for a particular simulation. Calculating the flux for either the fast 
cooling case 

Fu,max Vm>V>Vc (15) 

or slow cooling case 

i,^ / ^m) Pi/,max ^ ^ 

Fu={ {v/VraY^'^'^''''^ F^,max > 1^ > (16) 

i'^c/^my^^^^^^'^ {v/^cy^''^ Fu,max V > 

using the minimum Lorentz factor of the electrons 7^ = Ce (^fzf j comoving mag- 

netic field strength B = (327rmpeBn)^/^7c, and the cooling Lorentz factor of electrons 
7c = ^^B^^t — a^ifiR '^^^ ""^^ calculate the critical synchrotron frequencies and Vc 
using z/(7e) = 77e Standard values are taken for p = 2.4, eg = 0.1, = 0.01, n = 1 

(Panaitescu & Kumar 2002) and m^, me, cr-r and Qe are the proton mass, electron mass, 
fundamental charge and Thompson cross-section, respectively. Combining this with 

Fymax = — --fnBR^——, 17 

where D is the distance from the observer, here taken to be 10^^ cm, u is the frequency in the 
X-ray band, taken to be 10^^ Hz, one can calculate the X-ray flux density, F^, as a function 
of radius taking care to switch between fast and slow cooling cases where appropriate (Sari, 
Piran & Narayan 1998). 

Shells colliding onto the blast wave will produce "glitches" in the regular afterglow flux 
decay. Successive solutions, shown as dotted lines in Figl2], depend upon the energy contained 
in the blast wave and therefore on both the Lorentz factor, 70, and Mq, the effective mass 
of the blast wave. 



We use the calculated 7-evolution to calculate the evolution of the X-ray flux. In 
order to plot the flux as a function of time instead of radius, points are plotted as abscissa 
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t = i?/(2c7^) and ordinate flux where both flux and gamma are calculated by matching 
radius. Since the blast wave is accelerating in the rising part of Fig|2l here points appear 
to move back in time as the next solution is taken. In this part of the plot, blast wave 
solutions have simply been connected vertically. In the falling part of FigJ21 solutions appear 
to move forward in time as collisions occur. In general, all these abrupt jumps are artificial. 
In reality, one needs to consider the effect of equal arrival times, which smear out all the 
abrupt jumps so that the lightcurve would appear smoothed without noticeable individual 
glitches. 

Also shown in Figl2] are another blast wave solution for eg = 10"'^, and the afterglow 
template adopted in other calculations throughout the paper. As evident from the figure, 
for the standard value eg = 0.1 as derived from broad band afterglow modeling (e.g. Wijers 
& Galama 1999; Panaitescu & Kumar 2002; Yost et al. 2003), the X-ray flux predicted by 
the external shock model out-shines the template flux level (which is based on observations) 
by about three orders of magnitude. 

This is also emphasized in FigJSl which compares the simulated lightcurve (including 
the contributions from the prompt emission and X-ray flares superimposed on the template) 
and the calculated external shock lightcurve. It is evident that no steep decay phase and 
X-ray flares are observable if eg = 0.1 is adopted to calculate the afterglow level. This 
issue is carried over from the low efficiency problem of the internal shock model. If the 
observations are to be reproduced, there are two possibilities. The first is that the blast 
wave radiation efficiency is much lower. We test this possibility by lowering eg, and found 
that the predicted afterglow level can be roughly reproduced if eg is as low as 10~^ (Figsl2] 
and [3]). The pile-up effect of shells onto the blast wave naturally produces a shallow decay 
phase (after smoothing the abrupt jump features), which is an attractive feature of this 
internal-external-shock model. However, the anomalously small eg is inconsistent with the 
values derived from the previous broadband afterglow modeling, suggesting that this is likely 
not the correct approach to solve the problem. The second possibility is that the radiative 
efficiency of the prompt gamma-ray emission is very high (e.g. Zhang et al. 2007). This 
requires a more efficient mechanism to generate the prompt emission. The data of the recent 
GRB 080916C (Abdo et al. 2009) suggests that the outflow is very likely Poynting flux 
dominated (Zhang & Pe'er 2009). Within such a picture, Zhang & Yan (2009) proposed 
an Internal Collision-induced MAgnetic Reconnection and Turbulence (ICMART) model, 
which retains the merits of the internal shock model but significantly increases the prompt 
emission efficiency. 

The time history analysis of our shell model is also applicable to the ICMART model. 
The prompt emission and X-ray flare features can be retained, but the radiation efficiency 
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is increased. Within such a scenario, the external shock level can be lowered to satisfy the 
observational constraint. 

One can estimate the amplitude of "glitches" in the X-ray light curves. In general, 
collisions onto the blast wave are not energetic enough to produce a prominent signature 
on the afterglow light curve, unless the injection energy is comparable to that already in 
the blast wave (Zhang & Meszaros, 2002). In any case, small glitches from arriving shells, 
although not individually seen, could effectively bump up the normal decay phase, making 
it appear shallower. Postcollisions of shells sorted by decreasing Lorentz factor by internal 
collisions could produce an observable signature in early afterglow lightcurves and may be 
seen as temporal variability or a deviation from power-law decay (Kumar & Piran, 2000). The 
amplitude of ghtches produced by postcollisions can be calculated knowing the initial mass 
and Lorentz factor of both the blast wave and colliding shell and will essentially represent how 
far in the vertical direction successive solutions are separated. For the first few postcollisions, 
the amplitudes of the glitches are rather large as the added shells have both Lorentz factor 
and mass comparable to that of the blast wave itself. In this first stage, there is no simple 
approximation for calculating the jump in solution. The percent increases {{Ff — Fi)/Fi) 
are 3000%, 1000%, 500% and 150% for the presented simulation. During the deceleration 
phase, the X-ray flux density F^{X) oc ^(^+2)/^ for z/^ > min(z/^, z/J and F^{X) cx E^p+^'^I^ 
for Vm < i^x < i^ci where E is the total energy in the blast wave. The glitch amplitude is 
simply determined by the increase of the blast wave energy during each collision. During 
the collision, there is a reverse shock propagating to the trailing shell. However, this reverse 
shock usually does not contribute significantly into the X-ray band, since its density is higher, 
and hence, the typical electron Lorentz factor is much lower. Its dominant output is in the 
optical band (Meszaros & Rees 1997; Sari & Piran 1999; Zhang & Meszaros 2002). 

It has been argued that for GRBs with large peaks followed by deep troughs could be 
used to put an upper limit on the value of 70 (Zou & Piran, 2009). The external shock 
component will be superimposed on the internal shock prompt pulses and if this external 
shock does not rise above the level of the trough and/or threshold sensitivity of the detector 
at that point, this could be used to put an upper limit on the Lorentz factor of this first 
blast wave shell by knowing how high successive blast wave solutions will rise in flux. We 
find that this method, in principle, is a feasible method of putting a upper limit on 70 with 
a few caveats in mind. Figure [3] shows the blast wave solutions and blast wave (external 
shock component) plotted with the internal shock component lightcurve for one simulation. 
In general, the first shell of the external shock component will not be energetic enough to 
rise to the level of the prompt emission. As more and more shells pile onto the blast wave, 
the flux level rises and so looking for a solution that could be extrapolated backwards to find 
the initial 70 is problematic since the external shock component, especially early on, will not 
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0.1 10 1000 10^ 10^ 

t[s] 

Fig. 2. — The flux lightcurves are shown for the same blast wave used in Fig{T], calculated 
based on the external shock model. The calculation of eg = 0.1 is displayed in detail: the 
light dashed lines are the successive blast wave solutions, and the upper solid line is the 
corresponding blast wave lightcurve, which jumps between solutions. These glitches occur 
as shells pile onto the blast wave and decrease in magnitude as the energy of added shells 
becomes less signiflcant. The lower solid line is the lightcurve for eg = 10^'^, which matches 
the template adopted in the rest of the calculations (thick dashed line). 
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1 100 10^ 10^ 

t[s] 

Fig. 3. — The calculated X-ray lightcurve including prompt emission, X-ray flares and a 
broken power law template (solid black line), as compared with the lightcurves calculated 
from the external shock blast wave model (two colored lines: eg = 0.1 for the upper solid 
one and = 10~^ for the lower dashed one). 
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be a smooth function of time. Another comphcation is that the flux level of the blast wave 
emission sensitively depends on the unknown e^. parameter, making the derived upper limit 
of 7o subject to large uncertainty. 

4. Multiple Shell Simulations 

4.1. Single Injection Episode 

We first model collisions of a group a shells that are injected in a single emission episode. 
This is relevant to GRBs that have prompt emission without distinct gaps between pulses. 
Allowing the central engine to eject multiple randomized shells which go on to collide, we 
track the information of individual shells. An example is shown in FigsJl^ and E^. In 
this simulation, 100 shells are ejected where shell initial thickness, Lorentz factor and mass 
are chosen from random distributions in log space: Lorentz factors: 50 < 7 < 500, mass: 
10^^ < m < 10^^, initial thickness: 10^° < Aq < 2 x 10^°, all in cgs units. Ejection times from 
the central engine are chosen from a linear random distribution < t^j < 100 seconds in the 
rest frame of the GRB central engine. When two shells collide, we let them merge, drop one 
shell, and adopt the merged shell parameters (m, 7, A) as the new values of the remaining 
shell. In order to keep track of future collisions, we also need to re-set the "effective" ejection 
time of this new shell, which is taken as tej,m = Rcoi/cPm, where /3m = (1 — l/j'^Y^'^. The 
code then runs again with one shell reduced. The same procedure is applied when each 
collision happens, so that the code can track all the collision/merging processes for any 
arbitrarily designed central engine activity. 

Figure displays the "tree-plot" of this simulation. The upper panel displays how 
shells with different Lorentz factors are ejected, collide and merge at various distances. The 
time information is not displayed, but the collision times tcoi (again in the rest frame of 
the central engine) can be read off from the lower panel. This time is translated into the 
observed time according to Eq.Q to calculate the lightcurve (Fig|5^). The evolution of the 
blast wave is also marked in the upper panel as a thick red line. Collisions are "disallowed" 
if the shells collide with the blast wave (i.e. they would have collided at a radius greater 
than that of the blast wave, had the blast wave not existed.). These shells are included to 
boost the blast wave energy and become part of the blast wave after the collision and are 
therefore no longer traced in the later simulation. Our results suggest that a single episode 
can reproduce prompt emission of some GRBs. Even for simulations with a large number of 
shells, the ejection time in the frame of the central engine (tej) is correlated to the observed 
time of the collision (t©,co«) with little scatter (Fig. El see also Kobayashi et al. 1997). There 
are indeed collisions with much larger collision times (red crosses), but they are "excluded" 
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Fig. 4. — Tree plots showing simulated collisions as both Lorentz factor versus radius (top 
panes) and time of collisions versus radius of collisions (bottom panes). Lorentz factors are 
read up the y-axis starting at zero, time of collision is read down the y-axis, also starting 
at zero. The top pane shows shells as black lines, collisions are marked with dots. The 
time of collision can be found by dropping down vertically to the bottom pane of the graph, 
matching radius. The blast wave is shown in the top panels as a thick red line. Allowed 
collisions are blue dots, excluded collisions are red dots. For the left pane, 100 shells are 
ejected between and 100 seconds with Lorentz factors between 50 — 500. For the right 
pane, 50 A shells are ejected between < t^j < 100 seconds, 5 B shells 3000 < t^j < 5000 
and 5 C shells 30000 < ten < 50000 seconds. 
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by the blast wave constraint. If prompt emission and X-ray flares are indeed produced by 
the same mechanism, late flares are unlikely to be produced by shells that are ejected early. 
As discussed in the following, they demand re-activation of the central engine. Lightcurves 
are produced by considering the spectral model and temporal model as discussed in §2. The 
XRT band lightcurve of this particular simulation is shown in FiglS^. The relevant detector 
threshold (thick, dotted hne) is also added. (D. N. Burrows, 2008, private communication). 
It is evident that without late central engine activity, no X-ray flares can be detected. 

4.2. Multiple Injection Episodes 

In order to reproduce the observed X-ray flares that occur in distinct emission episodes 
from the prompt emission, we are obliged to simulate multiple ejection episodes from the 
central engine. The results are displayed in Figs. [71 Hb, [5)d and [H In these simulations, 
three groups of shells ("^", "i?" and "C") are released. The A group consisting of 50 shells 
is released in the flrst 100 seconds. The majority of these shells merge with each other 
producing AA type collisions. The energy released in each collision is shown as "AA" in 
Figd Similar to the previous simulation (§4.1), most A A type collisions occur before 100 
seconds. After 3000 seconds, we release a set of 5 i? group shells with energies reduced by 
~ 20 times on average. These shells can have two types of collisions: BB type collisions or 
AB type collisions (see FigjT]). Both the AB and BB type collisions could produce X-ray 
flares, as long as they are bright enough to stick out above the power law afterglow level 
(the broken power law template). For late X-ray flares, the energy of these collisions is of 
the "goldilocks type" : high enough to be above the background decay, but small enough to 
peak in the X-ray band and therefore remain undetectable by BAT. After 30,000 seconds, 
a third batch of C shells are released with energy lower than that of the B shells by ~ 20 
times again. This can in principle produce CC,AC or BC collisions, but for this particular 
simulation, AC and BC collisions are lacking (FigJT]). These cross collisions are disfavored 
since their collision radii tend to be large (due to the large ejection gap between the shells), 
so that the leading shells likely have collided onto the blast wave before the trailing shell 
catches up. Progressively less energy in A, B and C shells allows progressively degrading 
energetics of the X-ray flares, as is commonly observed. This is also generally consistent with 
various central engine models, where the accretion or magnetic power of the engine tends to 
die off with time. 

Another interesting topic is to investigate how the observed widths of X-ray flares can 
be reproduced. In particular, the observed 6t/t ~ 0.1 trend suggests that the later flares 
(larger t) are wider (larger 6t). This requires that the shell width A broadens with time. A 
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natural broadening mechanism is shell spreading (Eq.[7]). After a shell enters the spreading 
regime, the width of the shell is proportional to the radius, so that resultant X-ray flare 
width can be wide if the collision radius is large. Without placing the blast wave constraint, 
one indeed expects many "fat" flares, corresponding to very large radius collisions. With the 
blast wave constraint, the number of "fat" flares reduces significantly. Figure [9^ displays the 
model-predicted flare width as a function of their occurrences for a narrow distribution of the 
initial width Aq (between 10^° — 2 x 10^° cm, with the 6t/t = 0.1 line over plotted. The blue 
pluses are allowed but the red crosses are disallowed. The result suggests that although the 
predicted values are around the 6t/t = 0.1 line, the scatter is broader than what is observed. 
In particular, many narrow flares are predicted (those without significant spreading). The 
lightcurve for this simulation is presented in FigJTOk. which shows narrow late time flares 
that are not observed by the Swift/XRT data. In order to reproduce the observations, one is 
required to increase the initial shell width Aq for late ejection episodes. Figures (Hb andfTUb 
show an example of late injection of "fat" shells. Generally, Aq oc tej is needed to reproduce 
the observed data. This is consistent with the expectations of some central engine models. 
For example, in the fragmented disk model proposed by Perna et al. (2006), the clumps at 
larger radii have lower densities and tend to be more spread out so that the accretion time 
scale is longer. The ejected shells correspondingly also have longer durations. 



5. Conclusion 

We have developed a numerical code to model the internal collisions of an unsteady wind 
with arbitrary central engine activities. This is an extension of the internal shock models 
previously developed to model GRB prompt emission (e.g. Kobayashi et al. 1997), with a 
focus on X-ray flares that are commonly observed in GRB X-ray afterglows (e.g. Chincarini 
et al. 2007). Our motivation is to diagnose the required central engine activities based on 
the observational properties of X-ray flares as summarized in §1. The following conclusions 
can be reached: 

• The internal shock model with multiple ejection episodes can generally reproduce the 
properties of X-ray flares. Our shell model naturally explains both prompt emission 
and X-ray flares, suggesting that they originate from very similar mechanisms. 

• We find that the number of pulses/flares is directly related to the number of shells 
released from the central engine. More shells means more flares. In general, we find 
that the number of collisions is slightly smaller than the number of the shells ejected 
(e.g. the most probable values are: 5% less for 100 simulated shells, 10% less for 50 
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simulated shells, 8 collisions for 10 simulated shells and 1 collision for 4 simulated shells, 
see Fig|TT]for distributions). For a modal number of observed flares of 1, between 2 
and 5 shells need to be released on average. 

• The correlation between tej and te,co« ensures that the shells that are responsible for 
creating flares must be ejected just prior to being seen. This is because the second 
term in the right hand side of Eq.([H]) is typically much shorter than the first term. 
This means that not only the central engine activity must be prolonged, but it must 
also be episodic. Steady energy injection cannot produce flare-like features. In other 
words, early flares are created by shells ejected early, while late flares by shells ejected 
late. Since flares are seen as late as 10^ seconds, this means that the central engine can 
be active for a long time after the prompt emission. These conclusions are consistent 
with Zhang et al. (2006), Liang et al. (2006) and Lazzati & Perna (2007). 

• The large variance of fluences seen in flares can be explained by the different energies 
[E ~ jmc^) of the ejected shells. The peak luminosity of a flare also depends on the 
width of the pulse, which is either related to shell spreading or intrinsically different 
durations of shell ejection. 

• This study seems to rule out uniform (or narrow distribution) thickness shells. The 
observed widths of flares cannot be reproduced solely by spreading effects. One requires 
that later ejection episodes eject shells with larger widths to explain the typical flare 
width of 6te/t ~ 1/10. Since most shells collide before spreading, typical late collisions 
would be too narrow if later ejected shells had the same width as the initial prompt 
ejection episodes. Sharp flares where the width is near 1/10 of the emission time of the 
flare are common, but "fat" flares are seen occasionally as well. For example, the giant 
flare of GRB050502B has a 6tf,/t ~ 1 (Chincarini et al. 2007) implying that shells may 
either have spread before colliding or may simply have had a large width when ejected. 
Shells which have spread significantly before colliding will have energies spread over a 
large log- Gaussian shape and may be too dim to reach above the afterglow decay. 

• For goldilocks type shells, flares are seen in the XRT band only. The gamma-ray 
component of the X-ray flare (Band-function extension) is below the BAT detector 
threshold. Figure [8] shows light curves from a few typical simulations and these features 
can be seen. 

• Flares superimposed on other flares are simply shells that collide near the same time 
with different widths and fluences. The observed X-ray afterglow is a superposition of 
flares due to prolonged central engine activity and a background afterglow radiation, 
whose origin is not addressed in our paper. 
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• The decrease in average flare luminosity as a function of time indicates shells released 
at later times must create less energetic collisions (Lazzati et al. 2008). In order to 
produce the results seen, late released shells must be wider, less dense, slower or a 
combination of the above in order to create less energetic collisions seen. 

• The same shell model can give a prediction on the X-ray afterglow emission from the 
blast wave. Due to the continuous piling up of late-ejection shells onto the blast wave, 
the lightcurvc can show a shallow decay phase that is commonly observed. However, if 
a standard value = 0.1 is adopted, the predicted X-ray lightcurve is much brighter 
than what is seen, by about three orders of magnitudes. In order to reproduce the 
observed data, either a much lower (as low as 10~^) needs to be introduced, or the 
internal emission that powers the prompt emission and X-ray flares has to be much 
more efficient than internal shock predictions. 

Recently there has been interest in optical flares (Kriihler et al. 2009, Melandri et al. 
2006). These flares can be naturally interpreted in our model by invoking collisions between 
low energy shells or wide shells. Such collisions could be seen in optical bands but may be 
missed in higher energy bands. A direct expectation from this model is that optical flares 
should on average have lower energies and broader profiles than X-ray flares. An internal 
shock origin of optical flares was also proposed by Wei (2007). 

Finally, we want to emphasize that observationally flares have been seen in both Type 
I (e.g. GRB 050724) and type II GRBs. This requires that both types of progenitor have a 
similar central engine, which can eject an episodic wind to power late central engine activities 
which is a requirement of any central engine models for X-ray flares. 

We thank Shiho Kobayashi and Dave Burrows for discussion, and the anonymous ref- 
eree for helpful comments. This work is supported by NASA under grants NNX08AN24G, 
NNX08AE57A, AR9-0006X, and by NSF under grant AST-0908362. 
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Fig. 5. — A simulation with many shells in a single short episode will produce prompt 
emission only and no late collisions which could be responsible for X-ray flares. Left panel 
is a single ejection episode N = 100, < tej < 100 seconds. Right panel a simulation of N 
= 50 A shells, < tej < 100 seconds, 5 B shells 3 x 10^ < t^j < 5 x 10^ seconds and 5 C 
shells 3 X 10^ < tej < 5 X 10"^ seconds. Dashed lines show the position of flares underneath 
the superimposed afterglow. 
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Fig. 6. — The ejection time of a shell in the GRB central engine rest frame (tej) versus the 
collision time in the observer frame t©,co«- Allowed collisions are pluses, crosses are excluded. 
Even for a large number of shells (N = 100, ejected in 100 seconds), there is little scatter 
around a line of unity. 
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Fig. 7. — A bolometric light curve for 50 A shells ejected within the first 100 seconds, 
followed by 5 B shells after 3000 seconds and 5 C shells after 30,000 seconds representing 
a restarting of the central engine. Collision energies are shown by type; "AA" represents 
collisions between to A shells, "AB" between an A shell and a B shell etc. 



- 30 - 



Flux Light Curve in BAT band 



Flux Light Curve in XRT band 




1 10 100 1000 lO"* 10-'^ 10' 

tis] 

Flux Light Curve in BAT band 




10 100 1000 10-* 10^ 10' 

t[s] 

Flux Light Curve in BAT band 





Flux Light Curve in XRT band 




Flux Light Curve in XRT band 



1 10-' 

10"'' 







~^ \ — ■ ' 
, / 
' / 
1 / 





lo-* 

t[s] 



Fig. 8. — BAT (left) and XRT (right) light curves for three typical simulations of 50 A, 5 
B shells and 5 C shells of decreasing energy. Long dashed lines are the detector thresholds 
and short dashed lines in XRT lightcurves represent the positions of the flares underneath 
the afterglow. 
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Fig. 9. — How fat are the flares compared to tlie time of their occurrence? By the time 
most shells have spread, creating a broad flare, these shells have already collided with the 
blast wave or are not energetic enough to be seen. Observed flares are, on average, 1/10 of 
their occurrence times in width, this line is shown. Excluded collisions are shown as crosses, 
allowed collisions are pluses. The left panel is for shells of the same width, the right for shells 
which are wider for later ejection episodes. 
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Fig. 10. — Simulations with uniform shell widths (left) produce flares which are more narrow 
on average than observed flares. In order to reproduce observed wide flares, shells must be 
ejected with widths scaling as their ejection time (right panel). Notations are the same as 
for Fig. E 
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Fig. 11. — The most likely number of allowed collisions for simulations with 4, 10, 50 and 
100 shells. 1000 simulations were run for each. 



